/*
Este programa implementa o algoritmo de Borwein de convergência quártica, assim
como borwein.c. Este, porém, utiliza Pthreads para paralelizar parte do código;
especificamente, divide o trabalho em 2 threads, onde a segunda é responsável
pelo cálculo de b (independente de a e p a cada iteração), e a primeira é
responsável pelo restante. São utilizados dois semáforos: um antes do cálculo
de p na primeira thread (que é dependente de b) e outro para garantir que b
seja calculado somente após os valores da iteração anterior serem salvos.
*/

#include <pthread.h>
#include <semaphore.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <gmp.h>
#define NUM_THREADS 2

//Declaração das variáveis utilizadas nos cálculos
mpf_t a, b, p, pi;
mpf_t ap, bp, pp, aux[2][2];

//Semáforos e variável que indicará quando ocorrer a convergência
sem_t t1_loop, t2_loop, b_sem;
int end, i;

void* thread_1(void* dummy)
{
   while(1)
   {
      //Finalizar se convergir, e avisar thread 2
      if(mpf_cmp(p, pp) == 0)
      {
         end = 1;
         sem_post(&t2_loop);
         pthread_exit(NULL);
      }
      
      //Salva os valores atuais da sequência
      mpf_set(ap, a);
      mpf_set(bp, b);
      mpf_set(pp, p);
      
      //Sinaliza que thread 2 pode iniciar a iteração
      sem_post(&t2_loop);
      
      //Cálculo de a
      mpf_sqrt(aux[0][0], ap);
      mpf_ui_div(aux[0][1], 1, aux[0][0]);
      mpf_add(aux[0][1], aux[0][1], aux[0][0]);
      mpf_div_ui(a, aux[0][1], 2);
      
      //Espera pelo valor de b
      sem_wait(&t1_loop);
      
      //Cálculo de p
      mpf_add_ui(aux[0][0], a, 1);
      mpf_mul(aux[0][0], aux[0][0], pp);
      mpf_mul(aux[0][0], aux[0][0], b);
      mpf_add_ui(aux[0][1], b, 1);
      mpf_div(p, aux[0][0], aux[0][1]);
      
      //Incrementa o número de iterações
      i++;
   }
}

void* thread_2(void* dummy)
{
   while(1)
   {
      //Espera pela cópia dos valores da iteração atual
      sem_wait(&t2_loop);
      if(end)
         pthread_exit(NULL);
      
      //Cálculo de b
      mpf_sqrt(aux[1][0], ap);
      mpf_add_ui(aux[1][1], bp, 1);
      mpf_mul(aux[1][1], aux[1][1], aux[1][0]);
      mpf_add(aux[1][0], ap, bp);
      mpf_div(b, aux[1][1], aux[1][0]);
      
      //Sinaliza thread 1 que o cálculo de p pode ser realizado
      sem_post(&t1_loop);
   }        
}

int main(int argc, const char* argv[])
{
   //Declaração de variáveis para abertura de arquivo e medida de tempo
   FILE *fp;
   char filename[50];
   struct timeval ti, tf;
   long int sec, usec;
   
   //Início da contagem de tempo
   gettimeofday(&ti, NULL);
   
   //Declaração das threads e variáveis auxiliares
   pthread_t threads[NUM_THREADS];
   int t, precision, rc;
   void *status;
   
   //Inicialização de variáveis com a precisão especificada
   precision = atoi(argv[1]);
   mpf_set_default_prec(precision);
   mpf_init(a);
   mpf_init(b);
   mpf_init(p);
   mpf_init(pi);
   mpf_init(ap);
   mpf_init(bp);
   mpf_init(pp);
   mpf_init(aux[0][0]);
   mpf_init(aux[0][1]);
   mpf_init(aux[1][0]);
   mpf_init(aux[1][1]);   
   i = 0;
   
   //Valores iniciais da sequência
   mpf_sqrt_ui(a, 2);   
   mpf_add_ui(p, a, 2);
   
   //Inicialização dos semáforos
   sem_init(&t1_loop, 0, 0);
   sem_init(&t2_loop, 0, 0);
   end = 0;

   //Criação das threads
   rc = pthread_create(&threads[0], NULL, thread_1, NULL);
   if(rc)
   {
      printf("ERRO; codigo  da  pthread_create() e %d\n", rc);
      exit(-1);
   }
   rc = pthread_create(&threads[1], NULL, thread_2, NULL);
   if(rc)
   {
      printf("ERRO; codigo  da  pthread_create() e %d\n", rc);
      exit(-1);
   }
   
   //Aguarda as threads finalizarem
   for(t = 0; t < NUM_THREADS; t++)
      pthread_join(threads[t], &status);
   
   //p converge diretamente para pi
   mpf_set(pi, p);

   //Término da contagem de tempo
   gettimeofday(&tf, NULL);
   sec = tf.tv_sec - ti.tv_sec;
   if(tf.tv_usec < ti.tv_usec)
      sec--;
   usec = (tf.tv_usec - ti.tv_usec + 1000000) % 1000000;
   printf("Tempo decorrido: %ld.%06ld s. Iterações: %d\n", sec, usec, i);
   
   //Salva o resultado para arquivo
   sprintf(filename, "pi/borwein-par-%d", precision);
   fp = fopen(filename, "w");
   gmp_fprintf(fp, "%.Ff\n", pi);
   fclose(fp);

   return 0;
}
